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I derive the equation of motion in molecular dynamics for doing full lattice QCD 
simulations with clover quarks. The even-odd preconditioning technique, expected 
to significantly reduce the computational effort, is further developed for the simu- 
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1 Introduction 



The effects of dynamical quarks are important in QCD at finite temperature 
as well as in some phenomenological aspects at zero temperature. Unfortu- 
nately, the inclusion of dynamical quarks is the most demanding task in com- 
puter simulations of lattice QCD. The hybrid molecular dynamics or Hybrid 
Monte Carlo (HMC) methods have been developed into very efficient algo- 
rithms (maybe the most popular) for dynamical quarks. In these algorithms, 
the equation of motion is the essential ingredient. One has to derive the rel- 
evant equations before writing the programs for molecular dynamical simu- 
lations. For lattice QCD with staggered or Wilson fermions, these equations 
have already been available in the literature [1,2,3]. 

Lattice QCD has discretization errors due to the lattice spacing a. At inter- 
mediate bare coupling, corresponding to relatively large a, these systematic 
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errors might sometimes be very severe for Wilson fermions due to the chiral 
symmetry breaking term. The current computers do not allow the calculations 
done for very small a, because to reduce a implies to use a much larger lat- 
tice. Another way out is to use the improved fermionic actions. Recently, it 
has been shown that the use of the clover action [4,5] can significantly reduce 
these finite cut-off errors. However, the calculations of the clover action are 
much more complicated than the standard Wilson action. To my knowledge, 
there has not been a simulation of lattice QCD with the clover action in the 
literature. 

The purpose of this paper is to derive the equation of motion for full QCD sim- 
ulations with the clover action. Because the fermionic matrix has to be inverted 
in each step of the molecular dynamics step, it is also challenging to devise 
efficient algorithms for preconditioning [2,6] the fermionic matrix so that the 
inverse is easier to compute. For this reason, I also extend the even-odd pre- 
conditioning technique, previously used for quark propagator measurements, 
to the case of dynamical clover fermions. 



2 Preconditioning 

2.1 The action 

The action of the theory is S = S G + S F , where 



c p 



C X,fJ,>U 



£ Retrp^U^x + ^U^x + uYU^} 



(1) 



is the gauge action. The clover action for the quarks [4,5] is 
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B x,y = S(l - l^U^XjS^y^ + (1 + Jp)Ul(x - aO<Wm" (2) 

A is local and hermitian, and B connects only the nearest neighbor sites. The 
field strength tensor on the lattice is defined by T^ix) = [Q fll/ (x)—Q^ (x)]/2i, 
where Q^ v is the averaged sum of four plaquettes on the \w plane with the 
lattice site x as one corner. Each plaquette is the product of four link variables 
in the counterclockwise sense and begins with the link variable directed away 
from the site x and ends with the link variable directed towards site x, i.e., 



QA x ) = \[Up(x)u v (x + n)u^x + uyu u ( x y 



+U u (x)U tl (x - n + v) ] U v {x - n^U^x - //) 



+U tl (x - itfU v (x - At - v^U^x - At - v)U v (x - v) 



+ U u (x - u)%(x - v)U v {x + v)U^\ (3) 

as shown in Fig. 1. This operator is so chosen as for the maximal symmetry 
on the lattice. Most symbols in above equations are conventional, while the 
coefficient C in (2) depends on the choice of improvement strategy: C = 
1 for tree level improvement, and C = [Re tr(U p )/N c ]~ 3 ^ 4 for the tadpole 
improvement [7]. 



2.2 Even- odd splitting 



The lattice sites can be organized in an even-odd checkerboard and the even 
sites are numbered before the odd sites such that the fermionic matrix can be 
written as [8,9] 



M 



A ee kB co 



-kR a 



(4) 



/ 



where e or o denotes even or odd site on the lattice. 
Using such an arrangement, we obtain 
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Fig. 1. Plot of the clover operator Q^. The product in the plaquette is in counter- 
clockwise sense and begins with the directed link. 



det(M) = det(A 00 ) det (A 



ee 



n 2 B eo AjB oe ) = det A OQ det M ( 



ee 



(5) 



where 



M, 



ee 




(6) 



which couples only to even sites of the lattice. Now det(M) on the whole 
lattice has been factorized as a product of the determinant of the local matrix 
A on the odd lattice and that of M ee on the even lattice. 

To calculate the fermionic determinant, it is useful to introduce the pseudo- 
scalar variables r\ and <f) e , so that 




= J drjl drj expl-rjKAlAoo)- 1 ^} J d$ e d(j) e exp[-0i(M e t e M ee )- 1 0e], (7) 



where S p f is the pseudo-fermionic action 



S pf = rjKAlA^rjo + ^(j^m^-i^ 



(8) 



describing two flavor quarks with the same bare mass. Notice that r\ = A\ ^ 
and 4> e = M\ e d e have no direct coupling, where £ and 6 e are Gaussian noises 



injected at beginning of each molecular dynamics trajectory and held fixed 
during each trajectory. 

In the remaining text, I will use the above even-odd preconditioning to discuss 
the molecular dynamics. 

2.3 Fermionic inversion 

For quark propagator measurements and also in each molecular dynamics step, 
one has to calculate M" 1 or (Mj e M ee ) _1 , which can be implemented using the 
standard techniques like minimum residue, conjugate gradient or stabilized 
biconjugate gradient algorithms. The advantage of the even-odd splitting, as 
can also be seen later, is that such inversion is implemented only on the even 
lattice. Furthermore, due to the factor k 2 in (6), M ee is better conditioned 
than M. 

For the inversion A~l on each odd site, because it is completely local, we can 
use the LDV- decomposition [9,10] to solve it. Since it is a hermitian matrix, 
there exists a diagonal matrix D and lower-triangular matrix L such that 
A = LDL^ . Denoting % and j as the color-spin indexes of A, then 

i-1 

Di = An — E LikDkL* k , 
k=i 



i-i 

L H D j = A ij ~ E L ik D kL* jk i (.7 = 1,.. .,*-!)■ (9) 
k=i 

We can also compute the solution of AX = b by y = L~ x b and X = 
(L^D^Y, i.e., 

i-i 

Yi = bi - E L ikYk, (i = l,...,n), 
k=i 



X^Yi/Di- ]T L* ki X k , (i = n,..,l), (10) 

k=i+l 

with n = 12, which is the number of colors times the number of spins. The 
calculation is quite easy because there are n 2 multiplications and only n divi- 
sions. 
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3 Molecular dynamics 



3. 1 Equation of motion 



To develop the equation of motion for the gauge field A^x), one has to intro- 
duce a Hamiltonian H = J2 x ,fi Pa^(x)/^ + $g + ^p/> w ^ n ^M 1 ) tne canon i ca l 
conjugate momentum defined by P4 ( x ) = dH / d{dA IJt {x) / dr) , and r being the 
fictious molecular dynamics time. 

The gauge configurations are generated by solving the Hamiltonian equation 
of motion: 



d_Mx) _ dP A/l(x) dH dS G dS pf 

dr dr dApix) dA^x) dA^x)' 1 } 

For the gauge action, it is quite easy to show that 



dA^x) 2N c L dA fl (x) ^ ' MV ' dA^x) 



' ] [U^x)STAPLE^x)-h.c], (12) 



2N, 



where STAPLE^(x) is the sum over six staples surrounding the link U^(x). 
For the fermionic part, 



dS pf vA , d(Aj A 00 ) vA v , d{Mj e M ee ) v 
dA,(x) A ° dA,(x) A ° +Ae dA,(x) Ae ' [i6) 

where 



X* = (AlA^rjo, X e = (Mt e M ee )-Ve- (14) 
If we define two more variables 

Y A = A QO Xt Y e = M ee X e , (15) 
(13) can be simply rewritten as 
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+ Ae ^(x) ye + ye &4>) 



A straightforward computation leads to 



dM ee dA ee _ 2 ggeo 



i-re n eo /i OQ s± 00 n oe k n oe /i 00 



dMl dA ee dBj t 

cU M (z) dA^x) dA,(x) 00 oe 



+ K 2 Rt A- 1 dA °° A' 1 ^ -K 2 Rt 4" 1 dB ° e 
+ k B eo A OQ d ^ A OQ B oe k B oe A OQ 

By defining the following variables on the odd sites 

X Q = kA oq B oe X e , Y a = kA oq Bl e Y e , 
we have 



t dM ee _ f dA ee f <9A 00 



-^iw^-^iw> Y '- (19) 

We can further demonstrate that the last two terms of these two equations in 
(19) are summarized as 



KYe dA,(xf° KY °dA»{x) Xe KAe dA,(x) Y ° KA °dA,(x) Ye 



where 



= tr dirac [(l + !,)Y x+ ,Xl + (1 - ^)X X+ Xl (21) 

the same form as the fermionic force in the Wilson fermion case. As can also 
be seen later, the last two terms in (16) have exactly the same form for even 
and odd sites. Therefore, the introduction of the variables (14), (15) and (18) 
has a great advantage. 

A remark has to be made: to keep the conjugate momentum traceless, the 
right hand side (r.h.s.) of (11) should finally be subtracted by a term being 
the trace of the r.h.s. divided by N c . 



3.2 Practical implementation 



The simulations should be carried out in the first two steps as follows. 

1) Generating the full configurations. In numerical integration of the equa- 
tion of motion, one has to Taylor expand U(t + air) = exp[idrP(r)]U(T) = 
U{t) + idrP(T)U(r) + P{t + air) = P(r) + oIt^{t) + with finite 
order truncation. The leapfrog scheme can reduce the truncation errors to 
N m dO(dr 3 ) = 0(dr 2 ) at N m d molecular dynamical steps. These errors can be 
canceled by a Metropolis test at the end of the trajectory. Of course, this dr has 
to be fine tuned to maintain high acceptance rate and small auto-correlation 
time. 

2) About the clover coefficient. For the tree level improvement, C — 1. For 
the tadpole improvement scheme, the value of C depends dynamically on the 
configurations. One has to determine C self-consistently from the simulation. 
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For example, one may first have an initial guess for it, then generate a gauge 
configuration. From the plaquette, we can get a new C value. After some 
iteration, C might converge to some stable value for some given j3 and k. (This 
could be done very quickly since the plaquette can be accurately measured 
with a small number of configurations and on small lattices, provided the 
system is far away from a phase transition). 

3) Measuring the physical quantities. To obtain the improved hadronic matrix 
elements, rotation of quark fields [5,8] is necessary. 



4 More details about the fermionic force 



I have described that how the introduction of the variables (14), (15) and (18) 
leads to the equation of motion similar to that of the Wilson case. What is 
different is the terms with matrix A, which makes the equation of motion much 
more complicate. Therefore, it deserves further discussions. Note the term in 
the pseudo-fermion action Y^A 00 X^ (also the second term) is placed only on 
odd sites of the lattice, then for x being odd sites, there are terms only on x, 
x + /!+ v and x + /i — v relevant for Y^(dA 00 /dA^(x))X^ as shown in Fig. 
2, i.e., 



2 Y ° fe^dAM ° 



\E 'Y^a^[^^U u (x + f i)U,(x + uyU u (xy 



+U u (x - u)%(x - v)U v {x + /i - v)-^A^\xi 



1 - dU (xV 

+ 4 E ' y x %-^Mx + /i - v)-^L Uv {x - v)%(x - v)X^_ v . 
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x -v 



Fig. 2. Relevant plaquettes for Y] y [Yf{dA yy /dA^x))X^+X^{dAl y /dA^x))Y^\ 
when x e y, where the thick lines indicate the links relevant for dU^{x) / dA^{x) or 
dU^/dA^x). 



(22) 



Here J2l means the sum over v ^ fi. For x being even sites, there are terms 
only on x + //, x + v and x — v as shown in Fig. 3, i.e., 



' dA,{x) ° 



- ^ ' Y x A ^a, u [U v {x + ^(x + i/)tt^(x) 
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x +v • 



x O 



X -V 



Fig. 3. The same as Fig. 2 but for x doesn't belong to y. 



x+u 



+ \ E ' Y^Mx - v)U v {x + »- ")^^U„(x - v?Xt v 



(23) 



Therefore for odd-site x, the first two terms in (16) read 



yA\ 9A OQ A yAf ^oo yA 

° dA,(x) ° + ° dA^x) ° 
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o 



+ E ' F^^C/^x + uYU u (xyV^)U u (x + fi)X A 



- E ' Y x +»-v°Mx + fi - u)V^U v (x - v)%(x - u)X : 

v 

+ (Y ++X) + h.c.]. 
Similarly, for even x, the first two terms in (16) are 



yA\ 9A 00 A At dA\ Q A _ 
° dA,(x) ° + dA,(x) ° 



X + /1- 



+ E ' ^+U M ,C/,(^) t U M (x)C/ l/ (x + fiU^x + ^) f ^ +l/ 



E ' Y£ v <*Mx - u)U u {x + 11 - u)\J^U v {x - v^Xtv 



+ (Y ^ X) + h.c.}. 
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Equations (24) and (25) can be generalized to terms with y being odd or even: 



Viyt u - n yy y i yt vv vi 
Y y «dA»(x) Ay + A ydA»(x) Yyl 



= (24), if xey, 

= (25), i/a;^. (26) 
These relations are also quite useful when deriving the first two terms in (19). 
Summarizing (20) and (26), (16) becomes 

dS pf 
dA^x) 



= (25) + [(25) + (24)] (X A -> X, Y A -> F) + (20), i/ x = even, 



= (24) + [(25) + (24)] (X A -> X, F A -> F) + (20), i/ x = odd. (27) 

Note that the difference in the form of the molecular dynamics equation on 
even and odd sites is only in the first term. 



5 Discussions 

In this paper, I have derived (12) and (27), relevant for equation of motion (11) 
in molecular dynamics (or HMC) simulations with clover fermions. I have also 
extended the even-odd precondition technique, previously introduced for the 
quark propagator measurements, to the case of simulations with dynamical 
clover fermions. With the preconditioning technique, the number of iterations 
required would be reduced by a factor of 3 according to experience, and the 
most expensive part of the fermionic inversion is performed only on half lattice 
size. Therefore, it is expected that the preconditioned equation of motion 
would lead to considerable improvement over the unpreconditioned one. This 
scheme is vectorizable and has been parallelized. Currently, the simulations 
of QCD at finite temperature using the clover action are being carrying out 
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on the Quadrics-APEIOO parallel computers. The above work might lay a 
foundation of further computer simulations using dynamical clover fermions. 

Of course, to obtain physical results from the numerical simulations, there are 
still of lot of work to do. For instance, because the clover constant C depends 
dynamically on the gauge configuration, then there would be delicate interplay 
between C and [3, k, in particular when the system is at criticality. This is a 
new subject beyond the scope of the paper. 

It has been mentioned in [8] that even for the quenched clover propagator 
calculations, each minimum residue iteration took 35% longer than for the 
Wilson action. One has to choose a more efficient algorithm for fermionic 
inversion because a good algorithm is critical for full simulations. It is known 
that the stabilized biconjugate gradient is more efficient than the minimum 
residue or conjugate gradient, reducing the CPU time by at least a factor of 
1.5. 

Even if with the clover action there is an improvement of the finite cut-off 
error, the simulations with this action require larger statistical samples, more 
arithmetic operations and much memory. Concerning the feasibility of a full 
QCD clover simulation on supercomputers, it is not easy to quantify, because 
it is machine and code dependent. 

I hope to discuss these problems and report the physical results in the near 
future. 
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